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Abstract 

We present a novel way of accelerating hybrid surrogate methods for the calcula¬ 
tion of failure probabilities. The main idea is to use mesh refinement in order to 
obtain improved local surrogates of low computation cost to simulate on. These 
improved surrogates can reduce significantly the required number of evaluations 
of the exact model (which is the usual bottleneck of failure probability calcula¬ 
tions) . Meanwhile the effort on evaluations of surrogates is dramatically reduced 
by utilizing low order local surrogates. Numerical results of the application of 
the proposed approach in several examples of increasing complexity show the 
robustness, versatility and gain in efficiency of the method. 

Keywords: Failure probability, Adaptive mesh refinement, Multi-element, 
gPC. 


1. Introduction 

Most physical systems are inevitably affected by uncertainties due to natu¬ 
ral variabilities or incomplete knowledge about their governing laws. Over the 
past decade, the engineering research community has realized the importance of 
advanced stochastic simulation methods for reliability analysis. A quantity of 
paramount importance in reliability analysis is the calculation of failure proba¬ 
bilities. Mathematically speaking, it is a problem of computing multi-manifold 
integrals over some failure domains, whose structures are defined by some fail¬ 
ure functions, also know as failure modes. The irregular geometry of the failure 
domain, whose explicit form is usually unavailable, makes the accurate estima¬ 
tion of failure probability very difficult. In applications, the most interesting 
and relevant cases are usually in high dimensions and with very small failure 
probability, consequently the problem becomes more challenging. 

A straightforward method for estimating the failure probability (usually de¬ 
noted by Pf) is through Monte Carlo simulations (MC) (cf. [1]), where a num¬ 
ber of samples are drawn and the failure probability is estimated by counting 
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the proportion of the samples lying in the failure domain. The evaluation of 
the failure function (limit state function) is rather involved and the required 
sample size used in the evaluations is proportional to l/P/. As a result, the 
computational cost of brute force MC is usually prohibitively large. To alle¬ 
viate this inefficiency, various sampling techniques have been explored (e.g., 
low-discrepancy sampling [2, 3, 4], Latin hypercube sampling [5, 6, 7], im¬ 
portance sampling [8, 9, 10, 11], line sampling [12, 13], directional sampling 
[14, 15, 16, 17, 18], subset sampling [19, 20]). Other alternatives to reduce the 
simulation effort are non-sampling methods (which are also called approximated 
methods), such as FORM/SORM(first-order/second-order reliability method) 
(cf. [21, 22, 23, 24, 25]). FORM and SORM are based on Talyor expansions 
(first-order and second-order, respectively) of the failure function around the 
most probable point on standard Gaussian random space and use asymptotic 
analytic estimates for the failure probability. Although efficient, approximated 
methods have reduced accuracy because of the assorted approximating steps 
adopted in their implementations. 

The response surface method (RSM) is a combination of both simulation 
methods and approximated methods. To implement RSM, one constructs a 
surrogate (response surface) which is an approximation of the failure function 
from deterministic simulations on some design points or heuristic model, and 
then conducts MC on the obtained surrogate (cf. [26, 27, 28, 29, 30, 31, 32]). 
RSM has the advantage of the straightforward implementation of MC and that 
of the cheaper failure function evaluations on the response surface. However, as 
was shown in [33], RMS suffers from robustness issues. To address this problem, 
a hybrid method was proposed in [33], where both the response surface and the 
original underlying failure functions are utilized to construct the approximation 
of the failure probability. 

The method proposed in this paper combines this hybrid surrogate method 
with a recently developed mesh refinement method for random space [34, 35]. 
It is shown in [33] that a better surrogate results in more accurate estimates 
and requires less computation resources. As a result, to accelerate the hybrid 
surrogate implementation, constructing a better surrogate should lead to a more 
efficient algorithm. 

There are two ways that one can use to refine the resolution. One is p- 
refinement which amounts to increasing the degree of the polynomial basis 
functions used. The second is h-refinement which amounts to decreasing the 
size of the elements used. An accurate high order polynomial surrogate results 
in a lot of computational effort to be allocated in the evaluation of high order 
polynomials. However, if we keep the order of the polynomial basis functions 
only moderately large, we can at the same time refine the random elements 
into pieces such that the resulting multi-element polynomial expansion has the 
same error as the high order global surrogate. In the multi-element low order 
surrogate method, each sample point only needs to be evaluated by a low order 
polynomial while the multi-element decomposition maintains the accuracy of 
the surrogate. The adaptive mesh refinement methods proposed in [36, 34, 35] 
dynamically decompose the random space according to the evolution of the un- 
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derlying system. The algorithms automatically allocate more elements near the 
dynamical important regions which allows to maintain an adequate resolution 
everywhere. 

With this in mind, we propose two hybrid approaches which fit in the multi¬ 
element framework. In the first one we straightforwardly implement the iterative 
hybrid surrogate approach from [33] and replace the global polynomial surrogate 
by a multi-element surrogate. The second approach is recommended if one 
requires very high accuracy of the failure probability estimate. It implements 
the hybrid surrogate method in each element and then take the sum of the 
failure probability in each of the elements as the estimate. As the numerical 
results suggest, the second approach requires more evaluations of the exact 
failure function but also provides more accurate estimates. 

The rest of the paper is organized as follows. After presenting the formula¬ 
tion of failure probability computation in Section 2, we briefly review the key 
ingredients of the present method, including Monte Carlo simulation in Section 
3 and the hybrid surrogate approach [33] in Section 4. We rephrase the multi¬ 
element gPC expansion and adaptive mesh refinement methods as preliminaries 
in Section 5. The details of the new method are presented in Section 6, where 
we explain two versions of our approach. Results of numerical experiments are 
presented in Section 7 to demonstrate the effectiveness of our algorithms. 


2. Problem formulation 


Let (n. A, V) be a complete probability space, where fl is the event space and 
V is the probability measure defined on A G 2^, the a— algebra of subsets of fl 
(these subsets are called events). Let Z = (Zi, Z 2 ,..., ^ be a nz 

dimensional random vector with joint cumulative distribution function (c.d.f.) 
Fz{z) = Prob(Zi < ^ 1 , Z 2 < 2 : 2 , • • ■, Znz < z^), where 2 : = ( 2 : 1 , 2 : 2 ,..., z^) G 
A failure domain 12/ is usually defined by a limit state function g, also 
called performance function, failure mode or failure constraint as following: 

flf = {uj G fl : g(Z(oj)) < 0}. (1) 

Thus the failure probability can be defined as: 

Pf 4 Prob(a; G %) = f lz(Qp{z)dFziz), (2) 

dK"z 


where I is the characteristic function satisfying 


f 1 if z G A, 
\ 0 if z^A. 


( 3 ) 


Without loss of generality, we focus on the case ft C from now on, and let 
the joint probability density function(p.d.f.) of Z be q{z) = ■ 

The problem formulation is simple since the calculation of Pf amounts to 
the calculation of a nz—dimensional integral. However, Pj cannot be efficiently 


3 



calculated by direct numerical integration if nz is not small or the geometry of 
the failure region is not regular. The problem becomes more challenging if the 
failure probability is small. Since our method is essentially a mixed sampling 
method, we will briefly introduced the MC method. 

3. Monte Carlo Simulation 

The Monte Carlo simulation method is widely utilized due to its straight¬ 
forward implementation and robustness. Let = 1,... ,m be a set 

of samples of the random vector Z. The MC estimate of the failure probability 
is given by 

m 

i=l 

Note that, hereinafter we will, with a slight abuse of notation, use the short- 
handed notation {g(z) < 0} to stand for the set {z : g(z) < 0}. Although 
straightforwardly to implement, the MC approach can be costly in practice, 
since each sample point requires a full-scale simulation of the underlying sys¬ 
tem. Also, the convergence rate of this estimator is measured by the standard 
deviation of where = 1 /a/P/( 1 — P/)m [37]. Usually a large 
number of samples is required to obtain an accurate estimate of the failure 
probability, especially when Pf is small. 

4. Hybrid surrogate method 

A hybrid method was introduced in [33] to combine the robustness of MC 
with the efficiency of the surrogate approximation. We assume that there exists 
a surrogate model g that approximates the exact limit state function g. The 
surrogate can be constructed either by numerical approximations or through 
physicalmathematical reasoning. In general, we assume that the surrogate is an 
approximation of g with small norm error 

ep = ll5(Z) - ff(Z)|[LP = (J \g(z) - g{z)\Pq{z)dzy^P, p>l. (5) 

In [33] it was shown that straightforward sampling of the surrogate is fundamen¬ 
tally flawed, and may result in erroneous estimates, no matter how accurate the 
surrogate g is. To alleviate this disadvantage, a hybrid method was proposed, 
where the direct surrogate sampling result is corrected by results from the limit 
state function g. The key idea of the method is to replace the simulations of the 
surrogate g when the samples are ” close” to the response surface g = 0 by simu¬ 
lations of the limit state function g. This allows most of the sample simulations 
of the surrogate g to be kept. 
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The hybrid method seeks to estimate the failure probability using the fol¬ 
lowing two integrals, 


-P7 = y 

Q'f J' ^{—'y<g<'y}n{g<0}Q{z)dz^ 


( 6 ) 


where 7 > 0 is a small real number. With these integrals, the estimate of the 
failure probability by the hybrid method is 

P^=Py+Q^. ( 7 ) 

It was proved in [33] that with appropriately chosen threshold parameter 7 , 
the error of the hybrid estimate \P^ — Pf\ can be bounded by any prescribed 
accuracy control threshold e > 0. The choice of 7 depends on Cp, the norm 
{p> 1 ) of the error of the surrogate g, 

7 > ( 8 ) 

The hybrid estimate can be computed through MC as follows: Let be 

samples drawn from the distribution q{z), then 

1 ^ 

Pf = • (9) 

i=l 

Equation (9) requires knowledge of the prescribed threshold parameter 7 . For 
most practical applications the exact form of g is not available since one only 
has access to a black box model that provides a result for each sample. An 
iterative algorithm letting the algorithm automatically correct the estimate by 
gradually adding simulations from the exact model served as an alternative of 
the direct MC computation of the hybrid estimate. In the current work, we only 
implement the iterative version of the algorithm which avoids the prescription 
of the threshold parameter 7 . 


Iterative Hybrid Surrogate Algorithm 

Let TO 1 be the total number of samples, and S = 1 be the sample 

set generated according to the distribution q{z). Let Ato, called the ’’step 
size”, be an integer (much) smaller than to and 77 > 0 a small number used at 
the stopping criterion of the following iteration (also, let integer £ > 0 be the 
iteration count): 

• Initialization: 


— Set £ = 0, TO(^) = 0, and 57 ^) = 0. 
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~ Estimate the failure probability using the surrogate model g. That 
is, for i = 0, let 

m 

Pw = -T.h§<o}{z^% (10) 

Sort in ascending order. 

• Iteration: At t—th iteration (t > 0), do the following: 

— Identify the + 1 to (m(£) + Am) A m elements in the sorted 
sequence of and their corresponding samples in the set S. Denote 
AS'(f) the set of these samples, and let U AS'(^). 

— Evaluate the exact limit state function g at the sample points in 
A5(,). 

— Update the failure probability estimate using the value of g on AS'(f), 

P(i) = P(l-1) + — X] [“^{9<0}(^) +I{s<0}(^)]- (11) 

— If \P{t) — P(i-i)\ < r] or mg, + Am > m, exit; if not, let £•<—£+ 1, 
m(£) -(r- + Am, and repeat the iteration. 

• The failure probability is estimated by 

Pf = Pm- (12) 

In [33] it is shown that the accuracy of the surrogate directly impacts the effi¬ 
ciency and accuracy of the estimate of the iterative algorithm. To improve this 
hybrid algorithm, one approach is to improve the surrogate and another ap¬ 
proach is to improve the sampling strategy as presented in [11]. If the surrogate 
is constructed by polynomial expansion, we can use a higher order polynomial 
as the surrogate. Alternatively, we can divide the random space into elements 
and use a lower order multi-element polynomial expansion as the surrogate. 

5. Multi-element gPC surrogate 

For most problems of practical interest, the limit state function is obtain as 
the solution to an equation with random parameters, random initial conditions 
or random boundary conditions. We employ the generalized polynomial chaos 
(gPC) expansion to construct the surrogate for its exponential convergence rate 
when the system is smooth with respect to the random parameters. However, 
it is known that sometimes global polynomial expansion cannot capture the 
stochastic properties of this solution. In such cases, multi-element polynomial 
expansion solution will do a much better job. Since in hybrid surrogate method, 
we need the simulations of most samples on surrogate model, to gain the same 
accuracy low order multi-element polynomial surrogate is preferred to high order 
polynomial surrogate. In this section we will introduce the construction of multi¬ 
element polynomial surrogate via adaptive gPC mesh refinement. 
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5.1. Generalized polynomial chaos expansion 

The gPC method, an extension of the seminal work on polynomial chaos by 
R. Ghanem and (cf. [38]), has become one of the most widely used methods for 
stochastic computation since its introduction in [39]. Let i = • • • ,*nz) € 

Nq^ be a multi-index with |i|= -|- Z 2 H--I- inz, and let > 0 be an integer. 

The TVth-degree gPC expansion of g{Z) takes the following form 

N 

g^iZ) = Y, ( 13 ) 

|i|=0 


where {oi} are expansion coefficients that are to be determined, and {<i)i(Z)} 
are nz-dimensional orthonormal polynomials of degree up to N, satisfying 

E[$i(Z)ci>j(Z)] = [ $i(z)$j(z)g(z)dz = 0 < |i|, |j| < TV. (14) 

JQ. 

Here (5ij = Jldii ^id.jd = 1 if i = j aad 0 otherwise, and is the multivariate Kro- 
necker delta function. The orthogonality relation indicates a correspondence 
between the distribution of Z and the type of the orthogonal polynomial basis 
{<l>i(Z)}. For example, Hermite polynomial chaos corresponds to Gaussian dis¬ 
tribution, Legendre polynomial chaos corresponds to uniform distribution. For 
detailed discussions of these relations, see [39]. For a given stochastic system, 
there are different ways to obtain the expansion coefficients {ai}|)j^g such as 
stochastic Galerkin method and stochastic collocation method. We will not en¬ 
gage in a detailed discussion on this and refer the interested readers to references 
such as [40]. 

5.2. Multi-element decomposition in random space 

To deal with irregularities in random space, such as discontinuities, or to 
obtain a low order polynomial surrogate we discretize the random space Q C 
into a collection of non-overlapping hypercubes. Without loss of generality, let 
Z be a nz-dimensional random vector defined on the probability space H, where 
Zi,i = 1 ■■ ■ ,nz are independent identically distributed (i.i.d) random variables. 
Let fl be decomposed into M non-overlapping elements as follows: 

Bk = [oii&i) X [021^2) X ••• X [a^ 2 ’^^z)’ 

M 

n=\jBk, (15) 

k=l 

BiCi Bj = % if i ^ j, 
where i, j, = 1,2, • • • , M. 

For each random element the local random vector Z’^ = Z\\ has the 
conditional p.d.f 
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where = Prob(lBfc = 1 ) > 0 . 

Then a multi-element decomposition of the limit state function can be ex¬ 
pressed as: 

M 

g(Z) = ^5(Z)lB.(Z). (16) 

k=l 

5.3. Adaptive gPC mesh refinement 

Adaptive mesh refinement is a tool to dynamically refine the random space 
into multi-element according to the evolution or properties of the underlying 
system. From now on in this section, we solve the system in one element {B^) 
without specially indicating. Assume that the A^-th order gPC approximation 
of g can be stated as: 

N 

g{Z) = ^ 5 i<l>i(Z). (17) 

| i |-0 

Since ^ orthonormal basis with respect to the distribution of Z, the 

local variance is obtained by 


a 


2 _ 
N — 


N 


iii=i 


(18) 


If the system is steadied, mesh refinement criterion stated in [36] can be 
utilized to guide the refinement in random space. In this method, a local decay 
rate of relative error of the approximation in each random element is defined as: 


\i\=N 


(19) 


and the sensitivity of each random dimension is defined as: 


\i\=N 


( 20 ) 


where is the multi-index such that = 6ij, then Ne^ is the multi-index 
with j-th component N. Random element will be split when 

77 “Prob(Z G B'^) > 6»i, 

for some a G (0,1) and error control > 0. If nz = 1, will be divided into 
two equal parts, while if nz > 1 all dimensions that satisfy 

Ti > 02 max Tj, z = 1 ,..., nz, 
for some 0 < 6*2 < 1 will be split into two equal elements. 




If the system is dynamical and evolved with respect to time dynamical adap¬ 
tive mesh refinement algorithm stated in [34, 35] can be employed to get the 
multi-element gPC representation. Usually the dynamic of the system is ex¬ 
pressed as follows: 

ut{t,z) = C{t,z;u). (21) 

N 

Assume u = Ui{t)^i{z) is the gPC approximation of the solution to (21) in 
|i|=0 

one element. Then we obtain 

E E Mt)Mz)). (22) 

|i|=0 |i|=0 

Taking the Galerkin projection, (22) becomes 

( C{t,z]^u-^{t)^i{z))^i{z)q{z)dz, |i|= 0 , ...,Ai 

|j|=o (23) 

= -Ri(wj(0; |j| = 0,..., Af), 


here each .Ri(Mj; jjj = 0,..., N) denotes a function with uj, jjj = 0,..., as its 
variables, for jij = 0,..., N. If (23) has a reliable reduced model of resolvable 
variables Ui, jij = 0,..., A^o, Nq < N duS it is in [34]: 

^=i?i(uj(t);|j|=0,...,iVo), |i|=0,...,iVo. (24) 

In this case (23) is considered full system while (24) is the reduced system. We 
can compute the following quantity 


A^o o ~ Affl o ~ 


|j|=o 




|j|=o 

No 

= lE2i?i(wj(i);|j| = 0,...,iV)fii 

|i|=0 


No 

E2i?i(uj(t);|j|=0,...,iVo)ui|, 


|i|=0 


which can be considered as ’’the energy” transferred from the full system to the 
reduced system, (cf. [34]). The contribution of i-th dimension to Q is defined 
as: 


Si (Cij (t), |j I 0, . . . , 2Ii]sj^QZ (uj (t), |j I 0, . . . , Nq')U]sj^qz I. 

Random element will be split when 

QProb(Z e B'^) > 6»i, 

for some 6i > 0. If nz = 1, B^ will be divided into two equal parts, while if 
nz > 1 all dimensions that satisfy 

Si > 02 max Sj, * = !,...,nz, 
i=l,...,nz 
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for some 0 < 6*2 < 1 will be split into two equal elements. 

If the reduced system (24) is not available, as long as we have the Galerkin 
projection (23) of the system, we can consider the truncated system up to order 
Nq as the reduced model (cf. [35]) to do the adaptive mesh refinement. Ri, 

No 

|i| = 0 ,..., A^o is derived by plugging in u = in ( 21 ) and taking 

|i|=0 

the Galerkin projection. The mesh refinement algorithm is implemented when 
Q and Si are calculated with new set of 

After that we have a decomposition of random space and gPC approximation 
of the limit function (as function of solution to underlying system) in each ele¬ 
ments, then the multi-element gPC approximation of limit function is expressed 
as: 

M 

5(Z) = ^5(Z|b.)Ib.(Z) 

"m m n ( 25 ) 

« ^5'=(Z'=)Is.(Z) = E E 5f4>f(Z'^)Is.(Z). 

fe=l fc=l |i|=0 


6. Estimate of failure probability in multi-element framework 

Once we get the multi-element expression (25), the failure probability can 
be estimated as: 


P/ = / l!:ifqiz)dz = X! / 

fe=i 


(26) 


The quantity 

Pf = I Ia,nB^q{z)dz 

is the contribution to the failure probability from the element . 

Assume (17) is the approximation in element , then there exist two hybrid 
surrogate implementations. The first approach is to consider 

M 

g = ^g'=(Z'=)lB.(Z), (27) 

fe=i 

as the surrogate and implement the algorithm stated in Section 4. We call this 
approach multi-element surrogate global hybrid approach (ME-GHA). Alter¬ 
natively, since in each random element B^ we have an approximation 5 ^', the 
hybrid algorithm can be implemented in each B^ independently and provide an 
approximation . We call this approach multi-element surrogate local hybrid 
approach (ME-LHA). The ME-LHA algorithm follows 

Multi-element Surrogate Local Hybrid Iterative Algorithm 
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Let TO » 1 be the total number of samples, and S = be the sample 

set generated according to the distribution q{z). Let = S H C S he the 
samples that are located in the random element and for 

k = 1,... ,M, Let Am (called the ’’step size”) be an integer 

(much) smaller than m and 77 > 0 a small number used in the stopping criterion 
of the following iteration (also let integer i be the iteration count): 

• Loop over all elements fc = 1,..., M, 

— Initialization: 

+ Set ^ = 0, TO(^) = 0, and S{^t) = 0- 

* Estimate the failure probability contribution of the k-th element 
using the surrogate model . That is, for £ = 0, let 

= ( 28 ) 

i^l 

* Sort in ascending order. 

— Iteration: At th iteration (£ > 0), do the following: 

* Identify the TO(£) + 1 to (to^ + Am) A to^ elements in the sorted 

sequence of \g^\ and their corresponding samples in the set S^. 
Denote AS'(^) the set of these samples, and let S'(f) = U 

A5(,). 

* Evaluate the exact limit state function g at the sample points in 
+ Update the failure probability estimate using the value of g on 

AS'(£), 

-P(f) =-P(Vi) + ~ X! [“'^{s'“<0}(^) + “ll{s<0}(2)]- (29) 

zeAS(e) 

+ If — ^(^ 1)1 < 77 or TOf + Am > m^, 

Pf=Pw, (30) 

exit; if not, let £ ^ + 1, TO(f) ^ TO(£_i) + Am, and repeat the 
iteration. 

• End loop. 

• The failure probability is estimated as: 

M 

Pf=T.Pf- (31) 

k^l 
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Remark 1. There are multiple methods to construct the multi-element poly¬ 
nomial surrogates. We can simply divide the random space into elements and 
construct surrogates in each element as most of multi-element methods do no 
matter how the system evolves. Also, while we can adaptively refine the random 
space according to the evolution of the system based on different criteria. For 
example, one can use the contribution of highest order basis function to the total 
variance in one element [36]. Alternatively, one can utilize the rate of change of 
the activity transferred between full model and reduced model to decide on need 
for refinement [34] or the rate of change of energy between high order model 
and low order model in [35] and etc. In this paper, we use the latter approach 
to decide where to place the elements and how to construct the multi-element 
polynomial surrogates. 

Remark 2. Usually, the surrogates are constructed before implementing the 
failure probability algorithm. Compared to the sample-sized simulations, we don’t 
take the effort of constructing the surrogate model into account. As long as 
the surrogate is obtained, the computational cost of the hybrid method is mainly 
attributed to the simulations on exact model and evaluations of surrogates. Work 
for evaluation of the polynomial surrogate with the same order is the same. As 
a result, low order multi-element polynomial surrogates are preferable to high 
order global polynomials. This means that for the hybrid surrogate method, 
h—refinement is usually better than p—refinement. 

Remark 3. For ME-GHA, we directly implement the hybrid method on the 
multi-element surrogate, while for ME-LHA, the hybrid replacement of surrogate 
simulations by exact simulations takes place locally in each element. Since at 
least one iterative step is performed, if we take the same step size as ME—GHA, 
more exact sample simulations are required and thus, better accuracy is achieved. 
The step size is a prescribed integer, and we can choose a smaller one according 
to the problem to improve the efficiency. Also, if the total number of exact 
simulations is limited, an upper bound number of replacements can be set up as 
stated in [33[. 

7. Numerical Examples 

We present results of the hybrid surrogate approach under the framework 
of multi-element decomposition for a collection of examples of increasing diffi¬ 
culty: i) a simple example where a global surrogate will not work, ii) a simple 
linear ODE with a random parameter, iii) the Kraichnan-Orszag three-mode 
system with random initial conditions and iv) the Burgers equation with ran¬ 
dom boundary condition. 

7.1. A simple example 

In [33] the authors used counterexamples to demonstrate that direct Monte 
Carlo simulation of the surrogate model may result in erroneous estimate. One 
such example is when one considers a step function as the failure function and 
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global gPC expansions as the surrogates. In this subsection, we recalculate the 
failure probability for this example using a multi-element gPC surrogate. Let 
Z ^ [/[—1,1] be a random variable uniformly distributed on [—1,1]. Consider 
the failure function 

r -1, Z€[-1,0); 

g{Z) = { -0.5, Z = 0; (32) 

[ 0, Ze(o,i]. 

By construction g{Z) G L'^_i ip First, we consider the global gPC expansion 
for g using Legendre polynomials as basis functions. This approximation has 
the explicit form 


9^{Z) 


1 

2 



n=0 


(-l)”(4n-b3)(2n)! 

22"+2(n-hl)!(n)! 


P2n+l{Z), 


(33) 


where is the Legendre polynomial of order n (note that gP is the expansion 
of order 2p+ 1). Since {Pn}'^=o is a complete basis for L2_^ ^j, we have 

lim |lg-gP||i 2 =0. (34) 

p—>-oo 

It is well known that this expansion suffers from Gibbs oscillations. As a result, 
direct Monte Carlo simulation (MC) on the surrogates cannot give the correct 
estimate. A numerical simulation with 1,000,000 samples is implemented and 
the iterative step size in the hybrid algorithm is chose to be 1,000 due to the 
Gibbs phenomenon. The simulation of the same 1,000, 000 samples on the exact 
model is taken as the reference. 

Results of the direct simulation on surrogates gP{Z) of different p are pre¬ 
sented in the second row of Table 1. The results conhrm that in this case direct 
MC on surrogates show erroneous estimates even with very high order(15) global 
gPC expansion, while with the chosen step size, hybrid method gives the same 
results as the reference. Number of simulations on the exact model via hybrid 
method are shown in the third row of Table 1. Although hybrid method works 
in this case, half of the samples are have to be simulated by the exact model 
and thus, it is not efficient. 

Now consider the multi-element surrogate 

2 

9 = ^9^{Zi), (35) 

i=l 


where Zi = .^|[_i.o), ffi(-Z'i) = -1, for £ [-1,0), and Z 2 = 2'|[o4], 52 (-^ 2 ) = 
0, for Z 2 £ [0,1]. Since the multi-element polynomial expansion can exactly 
resolve the original limit state function, applying this surrogate, the hybrid 
algorithm converges after one iteration, which means only 1,000 simulations 
of exact model are required which is 1/1000 of the sample size compared to 
502/1000 of the sample size when applying hybrid algorithm on global gPC 
surrogates. The use of a multi-element surrogate for g can improve dramatically 
the efficiency of the hybrid implementation and at the same time maintain the 
accuracy of the estimate. 
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p 

0 

2 

7 

Prob(5P < 0) 

0.833187 

0.773777 

0.756490 

# 

502,000 

502,000 

502,000 


Table 1: Failure probability of step function estimate by global surrogates with different orders 


1.2. A linear one-dimensional ODE 

Consider the linear one-dimensional ODE 

^=-Zu, u{0]Z) = uo, (36) 

where Z is the random input. The failure function in this case is defined by 

f{u; Z) = u{T]Z) - Ud- 

The corresponding failure probability is defined as: 

Pf = Prob[/(M;Z) < 0], 

This equation has the exact solution 

u{t,Z) =uoe-^\ (37) 

Here we set all parameters the same as in [33] such that uq = 1, T = 1, Ud = 0.5 
and assume Z ~ is a Gaussian random variable with mean /i = —2 

and standard deviation <7 = 1. The reference is obtained by direct MC with 
1,000,000 samples. 

To construct the multi-element surrogate, we first express the random input 
as an expansion of a uiriform random variable X, X ^ [/(—1,1). Since $(a;) = 
i -I- ^erf)-^) is the cumulative distribution function (c.d.f.) of standard normal 
distribution, then 

^-I-V^crerf”^)^) ~ A/'(/r, ct). (38) 

Here erf(x) = Jo dt is the error function. Then Z can be approximated 
by an expansion of Legendre polynomials of a random variable X ~ [/(—1,1). 
We can write 

p 

Z^Y^hUiX), (39) 

i=0 

where Li is the Legendre polynomial of order i. 

We can apply the adaptive mesh refinement algorithm from [34] and we 
obtain a multi-element solution to (36). Since our surrogate is a polynomial 
of a uniform random variable, for comparison our MC samples are taken from 
uniform distribution and then transformed to a Gaussian random variable by 
(38) to get the reference probability Pref = 0.003541. The computational cost 
for simulation of a sample by multi-element polynomial surrogate of order p is 
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the same as the simulation of a global polynomial surrogate of order p. The 
advantage of our multi-element surrogate lies in employing only low order poly¬ 
nomials to obtain an adequately accurate surrogate. At the same time low order 
polynomials make the simulations of extremely large sample size affordable. 

Table 2 shows the results of the hybrid algorithm with multi-element surro¬ 
gate and the global surrogate. Given the same sequence of the sample points, 
each implementation resulted in the same outcome as the reference, however, 
the number of exact model simulations are different. The surrogate error has 



P 

3 

5 

7 

global surrogate 

# 

105,000 

54,700 

31,600 

number of elements 


5 

5 

4 

ME-GHA 

# 

3,700 

3,700 

900 

ME-LHA 

# 

4,100 

4,100 

1200 


Table 2: Results of two algorithm implemented on a linear ODE with multi-element polyno¬ 
mials surrogates of different orders p. 


mainly two components. One is the truncation error when we approximate the 
random input in (39). The second is the numerical error when we evolve the 
equation (36). Our adaptive mesh refinement method can reduce the second 
error. However to reduce the first error we need to increase p. To construct 
the multi-element surrogate, we fixed the tolerance {TOLi) which controls the 
quality of the mesh refinement. As a result, even with different order p, the 
value of the second error is similar for the two approaches. 

All implementations obtain the same failure probability estimate as the ref¬ 
erence. For the computational cost, when p = 3, the ratio of exact simulations 
in the multi-element implementation (ME-LHA) to that in the global surrogate 
implementation is about 1/26. The ratio is 1/13 when p = 5, while 1/26 when 
p = 7. Thus, we see that with the multi-element surrogates, the implemen¬ 
tation is much more efficient. As illustrated in the remarks, ME-LHA takes 
more samples of exact simulations than ME-GHA, and the more elements the 
larger difference of the number of exact simulations. In this case, ME-LHA and 
ME-GHA have the same accuracy. 


7.3. The Kraichnan-Orszag three-mode system 

Consider the following system obtained by a linear transformation performed 
on the original Kraichnan-Orszag three-mode system, 

dyi 
-JT = 


dy2 

dt 

dys 

dt 


-2/22/3, 

2 I 2 

—2/i + 2/3, 


(40) 
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subject to initial conditions 

2 /i( 0 ) = yi(0;a;), 2 / 2 ( 0 ) = y 2 ( 0 ; w), 2 / 3 ( 0 ) = 2 / 3 ( 0 ; w). (41) 

The solution of the system has a bifurcation depending on the value of the 
initial conditions 2/1 (0) and 2 / 2 ( 0 ). It was shown in [41] that a global Wiener- 
Hermite expansion cannot faithfully represent the dynamics of the system when 
the random inputs are Gaussian variables. Mesh refinement algorithms [36, 34] 
can efficiently quantify the uncertainty of the system when the random inputs 
are uniform random variables whose range includes the bifurcation point. In 
the current work, we examine only the case of one-dimensional random input. 

We choose the initial conditions 

2/i(0;w) = I, 2 / 2 ( 0 ; w) = 0.1,C(u;), 2/3(0;a;)=0, (42) 

where ^ ~ C/[—1,1]. The discontinuity (bifurcation) point 2/2 = 0 is contained 
in the random input space. We define the failure function as: 

f{uj) = yi{T;u;)-Ud. (43) 

The global gPC expansion cannot accurately describe the solution after about 
t = 8 . In our numerical experiments, we fixed T = 15 and Ud = 0.03. The 
reference result for the failure probability Pref equals 0.102651 which is obtained 
by MC with 1,000, 000 samples. The relative error is defined as: 

relative error = 

Pref 

First we employed the global gPC solutions in the hybrid approach with 
degree p = 3, 5, and 7. We found that the hybrid iterative algorithm cannot 
obtain a good estimate of the failure probability. In each implementation, the 
iteration stopped early due to the huge discrepancy between the exact solution 
and global gPC solution. As a result, the final estimates of the failure probability 
are not reliable with these global gPC surrogates. 

Next we tested the multi-element surrogates in the hybrid approach. By 
varying the value of the user-prescribed tolerance TOLi used to decide the 
need for refinement, we can obtain different multi-element surrogates. Table 3 
shows results for ME-GHA. Since the tolerance TOLi controls the error of the 
surrogate, the stricter the tolerance the better is the surrogate. This means that 
a smaller number of simulations from the exact model is needed when better 
accuracy of the estimates is achieved. 

If we fix TOLi, then by increasing p, usually we obtained a better surro¬ 
gate and fewer elements were needed. Table 4 shows results by ME-LHA. This 
means that we implement locally the hybrid method in each random element 
individually and sum up the contributions to the failure probability from each 
element to get the estimate of the failure probability. If we take the same step 
size as in ME-GHA, ME-LHA needs more simulations of the exact model than 
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ME-GHA because at least one iterative step is required in each element. How¬ 
ever these further steps increase the accuracy of the algorithm. Although in 
each implementation ME-LHA requires more exact simulations than ME-GHA, 
all implementations achieve the reference estimate. 



TOLi = lOe - 3 

TOLi = lOe - 4 

TOLi = lOe - 5 

p = 3 

No. of elements 

22 

38 

58 

# 

6900 

500 

200 

relative error 

0.06% 

0.16% 

0.015% 

p = 5 

No. of elements 

12 

22 

30 

# 

3900 

200 

200 

relative error 

0.012% 

0.14% 

0.021% 

p = 7 

No. of elements 

10 

16 

26 

# 

1400 

2200 

300 

relative error 

0.33% 

0.038% 

0 


Table 3: Results of ME-GHA for KO system approximated by polynomials with different 
orders p 



TOLi = lOe - 3 

TOLi = lOe - 4 

TOLi = lOe - 5 

p = 3 

No. of elements 

22 

38 

58 

# 

12245 

4700 

6029 

p = 5 

No. of elements 

12 

22 

30 

# 

3400 

3000 

3400 

p = 7 

No. of elements 

10 

16 

26 

# 

2900 

2200 

2800 


Table 4: Results of ME-LHA for KO system approximated by polynomials with different 
orders p 


7.4- Transition Layer for Burgers Equation 
Consider the viscous Burgers equation 

ut + uux = vuxx, 

u(—1) = 1 -I- (5, m(1) = —1. 


(44) 


Here 5 is a small perturbation of the left boundary {x = —1), and v is the 
viscosity. The solution to (44) has a transition layer, defined as the zero of the 
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solution at the steady state. The position of the transition layer is very sensitive 
with respect to the uncertainty 5 at the left boundary [42, 33]. Here we define 
our failure function as: 

f{z{5)) = -z{5) + zo, (45) 

where z(S) satisfies limt_>.oo u(t, z(S);S) = 0, and S ^ U(0, e), e <C 1. The value 
of z can be found as the solution to a nonlinear system of algebraic equations 

([42]). 

A A 

Atanh [—(1 + z)] = 1 + (5, Htanh[—(1 — z)] = 1. 

We use the criterion stated in [36] to adaptively construct the multi-element 
surrogate. Since we use a collocation method to construct the surrogate in each 
element, simulations of a new set of collocation points are required whenever 
refinement takes place (these simulations are counted as simulations of exact 
model in the cost of implementation). Here we use 21 collocation points in 
each element due to the high sensitivity of the transition layer with respect to 
the random input 5. We fix zq = 0.75 and calculate the failure probability 
Pf = Prob[/(z) < 0]. Reference failure probability Pref = 0.127478 is obtained 
by MC with 1,000,000 samples on exact model. We compared the results by 
global polynomial surrogate of different orders and corresponding multi-element 
polynomial surrogates in Table 5. While all the algorithms achieve the reference 
result, the efficiency is much better for the algorithm which uses low order 
surrogates (p=2,3). 

We compare ME-LHA to the global polynomial surrogate algorithm. We hnd 
that the ratio of the number of the exact simulations is 1/40 when p = 2, and 
1 /54 when p = 3. With global polynomial surrogate of order 5, 3921 simulations 
of exact model are required which is even larger than the number of exact 
simulations of multi-element surrogate of order 2. Evaluation on polynomials of 
order 2 is much cheaper than that on polynomials of order 5. 



P 

2 

3 

4 

5 

global polynomial 

# 

101,921 

62,921 

23,421 

3,921 

number of elements 


9 

7 

6 

5 

ME-GHA 

# 

1,757 

573 

431 

389 

ME-LHA 

# 

2,557 

1,173 

931 

799 


Table 5: Results of burgers equation via different surrogates (global polynomials and multi¬ 
element polynomials) with different orders 


8. Summary 

We have presented how the hybrid surrogate method for the calculation of 
failure probabilities from [33] can be modified to yield a multi-element random 
space decomposition. Two versions of the new algorithm were obtained: i) a 
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global hybrid algorithm with multi-element surrogates (ME-GHA) which is a 
direct implementation of original hybrid surrogate method [33] on multi-element 
surrogates and ii) a local hybrid algorithm of multi-element surrogates (ME- 
LHA) in which the hybrid method is implemented on each random element 
individually. The computational cost for both methods is much smaller than the 
original hybrid surrogate method. This is achieved firstly by the need for fewer 
exact model evaluations and secondly by the fact that each surrogate evaluation 
is cheaper since it uses lower order polynomials. Various numerical examples of 
increasing difficulty were used to demonstrate the efficacy and robustness of the 
multi-element hybrid surrogate algorithms. 
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